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Abstract. We introduce and analyze a virtual element method (VEM) for the Helmholtz problem 
with approximating spaces made of products of low order VEM functions and plane waves. We restrict 
ourselves to the 2D Helmholtz equation with impedance boundary conditions on the whole domain 
boundary. The main ingredients of the plane wave VEM scheme are: i) a low frequency space made of 
VEM functions, whose basis functions are not explicitly computed in the element interiors; ii) a proper 
local projection operator onto the high-frequency space, made of plane waves; in) an approximate 
stabilization term. A convergence result for the h -version of the method is proved, and numerical 
results testing its performance on general polygonal meshes are presented. 
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Introduction 

The Virtual Element Method (in short, VEM) is a generalization of the classical finite element method 
recently introduced in HE]. The key features of the VEM are: 

• the decomposition of the computational domain can consist of arbitrary polygons in 2D and polyhedra 
in 3D; 

• the local spaces, in general, contain polynomials, as in the classical finite elements, but include also 
more general functions (usually defined through a differential operator) whose pointwise values need 
not to be computed (hence the name “Virtual”); 

• the VEM passes the Patch Test. 

The basic VEM introduced on the Poisson problem has then been extended to highly regular m , non con¬ 
forming p] and discontinuous approximating spaces [16], to general elliptic [8] and mixed problems [BJ[7J[2] , as 
well as to Stokes pQ, plate problems EE linear and non linear elasticity [SJIIDIGS] , and fluid flows in fractured 
media [12]. 

In this paper, we aim at extending the VEM approach to an indefinite problem. More precisely, we present 
and analyze a new method, based on inserting plane wave basis functions within the VEM framework, in order 
to construct a conforming, high order method for the discretization of the Helmholtz equation. As in the 
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partition of unity method (PUM, see e.g., [3TH39] ). we use approximation spaces made by products of functions 
that constitute a partition of unity and plane waves. 

Plane wave functions are a particular case of Trefftz functions for the Helmholtz problem, i.e., functions 
belonging to the kernel of the Helmholtz operator. Other Helmholtz-Trefftz functions are, e.g., circular/spherical 
waves (also denoted as Fourier-Bessel functions, or generalized harmonic polynomials), or Hankel functions. The 
idea of inserting Trefftz basis functions within the approximating spaces in finite element discretizations of the 
Helmholtz problem is due to the fact that these spaces possess better approximation properties for Helmholtz 
solutions, compared to standard polynomial spaces (see, e.g., [42]), thus similar accuracy can be obtained with 
less degrees of freedom. This mitigates the strong requirements in terms of number of degrees of freedom per 
wavelength due to pollution effects [3] . 

There are in the literature several finite element methods for the Helmholtz problem which make use of Trefftz 
functions. Besides the already mentioned PUM, which is H 1 -conforming, other approaches use discontinuous 
Trefftz basis functions and impose interelement continuity with different strategies: by least square formulations 
(see [53] 06]); within a discontinuous Galerkin (DG) framework (like the ultra weak variational formulation 
(UWVF) p~8lH9] or its Trefftz-DG generalization [281I3QU32] : see also [171124] for the derivation of the UWVF as 
a Trefftz-DG method); by the use of Lagrange multipliers (see, e.g., [22]|23]|33]|3H] ); through weighted residual 
formulations (like in the variational theory of complex rays [35J05], or in the wave based method |201l21| h 

As a first construction of a plane wave-virtual element method (PW-VEM), we focus here on the 2D Helmholtz 
problem with impedance boundary conditions on the whole domain boundary. We restrict ourselves to low order 
VEM functions and we take a uniform plane wave enrichment of the approximating spaces. On quite general 
polygonal meshes, our basis functions will be then products of low order VEM functions associated to the mesh 
vertices, multiplied by a linear combination of p plane waves centered at the corresponding vertices. 

The VEM framework is not only used for the choice of the low frequency space , whose basis functions are 
not explicitly computed in the elements interiors, but also for the definition of the method and the consequent 
approach to its analysis. A crucial step is the definition of a proper local projection operator onto a space 
that has to verify two major requirements: providing good approximation properties for the solution of the 
homogeneous Helmholtz problem, and allowing to compute exactly the bilinear form whenever one of the two 
entries belongs to that space (consistency requirement). For this reason, we define the local projector through 
the Helmholtz bilinear form onto the space of discontinuous plane wave functions. 

The discrete bilinear form is then split into two parts: one that can be computed exactly, thanks to the 
projection operator, and a second one with the role of a stabilization term and that can be approximated. 

The analysis of the method is given in an abstract form and a sufficient condition on the stabilization term 
that implies the discrete Garding inequality (which in turn implies convergence) is provided. 

The outline of the paper is as follows. We introduce the model problem and its PW-VEM discretization in 
Section [I] describing the approximating spaces and the projection operator, as well as all the relevant matrix 
blocks needed for the implementation; the choice of the stabilization operator is not specified at this stage. 
In Section [2] we prove an abstract convergence result of the h -version of the PW-VEM which, together with 
best approximation estimates for the considered approximating spaces, give convergence rates of the method, 
provided that continuity and Garding inequality are satisfied for the discrete operator. A sufficient condition 
on the stabilization term that guarantees a Garding inequality for the discrete operator is given in Section [3] 
An explicit form of a possible stabilization term is then provided. Finally, we numerically test in Section 0] the 
performance of the resulting PW-VEM. 


1. The PW-VEM method 

We introduce the PW-VEM method for the Helmholtz problem. For simplicity, we consider the two- 
dimensional case. 
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Let H C M 2 be a bounded, convex polygon. Consider the Helmholtz boundary value problem 


—A u — k 2 u = 0 in (7, 
Wu • v + iku = g on dfl, 


(1) 


where the unknown u is a complex-valued function, k > 0 is a given wave number (the corresponding wavelength 
is A = 2ir/k), v is the outer normal unit vector to dfl, and i is the imaginary unit. To simplify the presentation, 
we are considering impedance boundary conditions on the whole dd, with datum g E L 2 (dfl). 

Since we are interested in the mid- and high-frequency regimes, we assume that A < diam(fl) or, equivalently, 
k > 2 tt/ diam(fi). 

The variational formulation of m reads as follows: find u E H 1 ^) such that 




( 2 ) 


where 



see, e.g., [37, Prop. 8.1.3] for the well-posedness of the variational problem (|2]). 


1.1. Discretization spaces 


Let Th = {K} be a mesh of the domain Q made of polygons K with mesh width h, i.e., h = ma xxeT h hx, 
where hx is the diameter of K. Given K E Th, we denote by the mass center of K , and by vx the unit 
normal vector to dK pointing outside K. Moreover, we denote by nx the number of edges e of K, and by Vj 
and Xj, j = 1,..., riK, the vertices of K and their coordinates, respectively. 


Let K be an element of Th- First of all we choose the low frequency space as the (local) finite dimensional 
space V(K) defined as 


V{K) = {ve H\K),v ldK E C°(dK), v\ e E Pi(e) Me C dK, Av = 0 in K}. 


(3) 


Here and in the following, we denote by Pi(D) the space of polynomials of degree at most one in the domain D. 
The functions in V(K) are completely determined by their value in the uk vertices. The dimension of V(K) is 
equal to We denote by the canonical basis functions in V(K) defined by 


— i,j = 


Notice that ¥i(K) C V(K), and ¥i(K) = V(K) if and only if uk = 3. 

The basis is a partition of unity, i.e., /(x) = = 1 f° r all x E if. In fact, since the 

functions ipj are harmonic in K , linear on each e C dK , with (pj(Vi ) = Sij , then A/ = 0 in K and / = 1 on 
dK , which imply / = 1 in K. 

Next, we introduce px different directions {d^}^, and we define the local PW-VEM space 


n K Pk 


V P k( k ) = \v: v = Y^Y^a j <np j (x)e zkdt ' {x x ^, a je G c}. 


v : v = 


j =i £=i 


We also introduce the standard plane wave space (centered at xx) 


PK 
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Clearly, V p * K (K) C V Pk (K). 

Setting, for t = 1,.. * and j = 1,..., 

p«;/(x) = e ifedf '( x - x *), pwjg (x) = e ikde < x ~ x i\ 
ipr(x) = (pj (x) pWjg (x) , with r = (j - 1 )p K +1, 


we have 

Ux (if) = span{Vv(x), r = 1,..., n^px}- 

Notice that pWjf (x) = c :l c pu!f(x), where Cjp = e lM <'( x K _x i). 

Finally, we define the global PW-VEM space 

V p (%) = {»£ C°(U) : v\k e V P (K) VK e %}, 

where we have chosen pk = p and the same directions {d^}^ =1 for all K E 7^,, which allows to impose continuity 
across interelement boundaries. 

By discretizing m in the spaces V p (7h) one obtain the PUM method [5T1I59] . The space V(K) is the 
polygonal finite element space with harmonic barycentric coordinates (see, e.g., m), that can be considered 
also as VEM space of lowest order; see m for other choices of generalized barycentric coordinates. Here, we 
adopt the VEM framework and avoid quadrature on polygons, as well as the expression of the basis functions 
(fj in the element interiors. 

For the validity of approximation estimates in plane wave spaces (see [411142] and the proof of ProDosition ll.il 
below), we make the following assumptions on the meshes and on the plane wave directions: 

i) there exist p E (0,1/2] and 0 < po < p such that every K E Th contains a ball of radius phx and is 
star-shaped with respect to a ball of radius poh; 
ii) there exists an integer m > 1 such that 

p = 2 m + 1; 

in) the directions satisfy a minimum angle condition, i.e., there exists 0 < S < 1 such that the 

minimum angle between two different directions is > ^5, and are such that the angle between two 
subsequent directions is < 7r. 

1.2. The projector II 

In the VEM framework, a crucial step is the choice of a local projection operator that allows to compute 
bilinear forms without the need of having at disposal the explicit form of the basis functions (the functions ijj r 
here). The space where projecting onto has to fulfill two major requirements: providing good approximation 
properties for the solution of the problem at hand, and allowing to compute exactly the bilinear form whenever 
one of the two entries belongs to that space. 

Here, we choose as space where projecting onto the space of (discontinuous) piecewise plane waves, which we 
denote by i.e., 

v;{T h )= n v;(K). 

KeT h 

For any K E Th, we define the local bilinear form 

a K (u,v)= j Vu-VvdV-k 2 f uvdV, 

Jk Jk 

which is a local version of the bilinear form defining the variational problem ([2]), ignoring the boundary term. 
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We define the projector II : V P (K) —>> V*(K) as follows: 

a K (Uu,w) = a K (u,w) Mw G V*(K). (4) 

The projector II is well-defined, provided that k 2 is not a Neumann-Laplace eigenvalue on K. This is guaranteed, 
for instance, when condition ([5]) below is satisfied (see Section ip for more details). 

Clearly, if u G V*(K), then Hu = u. We observe that the choice of the low frequency space V(K) and of 
the high frequency space V* (K) allows us to compute the right hand side of ([!]), even if we do not know the 
expression of the functions of V(K) in the interior of K. In fact, for any u G V P (K) and w G V*(K), integrating 
by parts, we have 

a K (u,w)= / Vr • \7w dV — k 2 / uwdV = / u\7w • vk dS, 

Jk Jk JdK 

since Aw + k 2 w = 0. The integral on dK can be computed because both u and Ww • vk are known on dK. 

1.3. Continuity of the operator II 

Let fi 2 be the smallest strictly positive eigenvalue of the Neumann-Laplace operator on K. If K is convex, 
then fi 2 > t 2 /h 2 K [44]. If K is star-shaped with respect to a ball, /i 2 > CoiY 2 /h 2 K , with 0 < Co < 1 only 
depending on the shape of K USE]. Therefore, assuming that 

0 < hxk < Ci < y/Co7r, (5) 


we guarantee that k 2 < /i 2 - 

We denote by ll'llo the L 2 -norm in the domain D , and define the weighted norm 


\\ v \\i,k,K = \\^ v \\o,k + k2 Mo,k VveH\K). 

The following discrete inf-sup condition holds true. 

Proposition 1.1. Provided that hxk satisfies 


hx;k < ao < min 


{ 


\/Co 7T 

T7T 


,0.5538}, 


(6) 


there exists a positive constant {3 = ^(hxk) (i.e., /3 depends on hx and k through their product hxk), which 
remains uniformly bounded away from zero as hxk —»• 0, such that 


VveV;{K) Bw G V* (K) s.t. 


a K (v , w) 


> f3(hxk) IMIi^x • 


(7) 


Proof. Assume, with no loss of generality, that the direction d* = (1, 0) belongs to the set of directions {d^}^ =1 . 
Due to our assumptions on the directions (see assumption Hi) at the end of Section [Tj) , there are two directions 
d^dg G {di}^ =1 such that d^d^dj are listed in counterclockwise order, and the angles 5 i, 62, £3 between dg 
and d^, d^, and dj, d^ and dg, respectively, satisfy 0 < sin$i, sin £ 2 , sin $3 < 1 . 

Let b\ G V*(K) be defined as 

3 

b\ (x) 

1=1 


with 


sin £3 

sin S 1 + sin S 2 + sin $3 ’ 


sin Si 

sin Ji + sin S 2 + sin $3 ’ 


sin S 2 

sin <5i + sin S 2 + sin $3 


a 1 — 


OL 2 = 


& 3 = 
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Since Yle=i a £ = 1? then fri( x if) = 1, and since Yle=i a £^i = then V6i(xx) = 0. 


The following bounds hold true: 


IMIo,*<l*1 1/2 , 


I|l-&i|| 0 ,«:< ^Kk)*\K\V 2 , 
< h K k 2 \K\^ 2 . 

The bound (J5| follows from a \, 0 : 2 . <'>'3 > 0 and YJe=i a t = 1 ; 


(8) 

(9) 

( 10 ) 


IIMiU = j K Ml' dV ^ j R (X>)V 


= \K\. 


The bounds © and (non immediately follow from Taylor’s expansions of &i(x) and V&i(x), taking into account 
that b\ (xx) = 1 and V6i(x^) = 0. 

Fix 0 / v G V* and define w as the unique element of V* such that 


L 


■/. 


Ww-V£dV + k 2 w£dV = a K (v,£) V£ G V* 


( 11 ) 


Ik Jk 

Since the threshold condition (j5j) implies that k 2 < the right-hand side is a nonzero functional of £, therefore 
w / 0. Notice that, taking £ = v — w in (jTTj) , we have \\\/{v — w )||q K = k 2 ||v||q k — k 2 ||ic||q k which implies 


IIV('L' w) 11^ IMIo,if • 

Taking b\ as test function in (fill) , we have 


( 12 ) 


L 


vbidV = - 

k Jk 

(fToli 


f wbidV + -^ ( {Vv -Vw) -VhdV 
Jk k Jk 


< ~ f wb 1 dV-^\\V{v-w)\\ 0K h K \K\ 

Jk 


1/2 


(13) 


(fl2l) 

< 


• [ whdV + kWvWnjchKlKl 
Jk 


1/2 


Moreover, by taking £ = w in m, we have 


a (v,w) = |Mli, fc ,if 

In order to conclude, we need to prove that || / u || 1 ^ ^ ||w|li,fc,if • 

From the definition of a x (-, •) and from cn, setting c v = ^ f K v dV , we have 

IMIi /C if = aK ( v ’ v ) + 2 ^ 2 Iloilo k = [ ViT • Vpdy + k 2 f wvdV+ 2k 2 \\v\\l K 

’ J k Jk 

— \\ W \\l,k,K W V Wl,k,K + ^ \\ v ~ c v\\o,K + ^ Iloilo, if • 

The min-max principle implies 


(14) 


I|V<k ^ 


T - c, 


NO,if 


6 if 


2 k \\v c v 11o,if — 


a /2 hxk 


l|Vi 


lo,if • 


(15) 
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For the term 2k 2 ||c/|q k , we have 

k 

W\ 175 

k 


k II c 


•^Il0 ,K 


(fl3l) 

< 


< 

(EJ,11) 

< 


L 

L 


vdV 


wh dV 


| Kl 1 / 2 


f v b\ dV + f u(l-fti) 
Jk Jk 


dV 


+ (h K k)k ||«Ho,*r + Tk^72 J K V(l-h)dV 


\K\V 2 

Iloilo, K ll^lllo.K + (hxk)k Iloilo K + I 7^11/2 Iloilo, K I|1 — Mlo./f 


\K\V 2 


\K\ 1/2 


(16) 


k Ikllo.x + (h K k)k |M| 0)Jf + |r|| (UV 


— ^ Iloilo ,K + 

Inserting (USD and m into m gives 


1 + 


hxk 


(h K k)k (ML 


K • 


IMIl,fc,K — \\ W \\l,k,K W V \\l,k,K ' 




hxk 


1 2 


(WU Iloilo,*: 


+4 


hxk 


( hxk)k ll'R’Ho k ^ IMIo,k 


— + 4(hx/c) + 2(hi^k) 2 ) ll'R’Hi i. k IMIi, 


k,K 


a/2 hxk 
v / Co" 7r 


1 + 


hxk 


n 2 


(^) 2 


where in the last step we also have used k 2 ||w||q k < k ||ie|| 0 K k ||t|| 0 k \ thus 


1 — max 


{(^) 2 ’d! + ¥f (w 2 } 

3 + 4 (tiKk) + 2 (JiKk) 2 


\l,k,K — IMIl,fc,K' 5 


with a positive coefficient on the left-hand side, provided that max | ( ) 5 2 [l + (hx&) 2 | 

’ is convex, for instance, Co = 1, and th 

- max { (^) 2 ’ 2 [! + ^] 2 (U^) 2 


< 1, 


i.e., provided that (|6j) is satisfied (if K is convex, for instance, Co = 1, and the condition is hi^k < 0.5538). 
Then, (J7]) holds with 


P{h K k) = 


3 + 4 (hxk) + 2 (hxk) 2 


□ 


Remark 1.2. In the infinite dimensional case, under the restriction 0 < hxk < a\ < v ^| 7r , the following 
inf-sup condition holds true: 

\/v e H\K) 3w e H\K ) s.t. > P*(h K k) \\v\\ 1XK , 


W 


l,k,K 


with (3*(h K k) = 1 - ■ This can be proved along the same lines as in the proof of Proposition E 

choosing b\ = 1 since, in this case, it is an admissible test function. 
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We have carried out numerical tests, in the convex case (Co = 1) in order to numerically determine the inf- 
sup constant in with only the above restriction on the product hxk. The results obtained seem to indicate 
that the function /3 is bounded from below by the function /3*(hKk), and thus that the stronger restriction (|6]) 
on hxh we have required in Provo sition \l.l\ might not be needed in practice. 

The following result is a straightforward consequence of the inf-sup condition ([7]). 

Proposition 1.3. Under the condition the operator II is well-defined, and the following local continuity 
continuity property holds true: 


< /3(hKk) ^ Vp(K)- 

1.4. The matrices D , T>, G, and P 

With the same notation as [9], we introduce the basic components of our PW-VEM. 

Let D be the matrix of size ( ukP , p) whose £-ih column contains the coefficients of the representation of the 
plane wave pwi in the V P (K) basis {Vv}- Simple calculations, taking into account that is a partition 

of unity, show that the only entries of D different from zero are 

D((j - 1 )p + 1,1) =* e- ifed H x -- x y j = 1 ,..., n K , £ = 1, • • • ,p. 

We define B as the matrix of dimension (p, ukp) such that 

B(£,r) = £ml ,... ,p, r = 1,... ,n K p- 

As already observed in the definition of II, since A pw^ + k 2 pw^ = 0, the entries of B can be computed in 
terms of the traces of the shape functions on dK only: 

B(£,r) = a K ('ipr^pwi) = / \7'ip r ■ \7pw £ dV — k 2 / 'iprpw^dV 
Jk Jk 

= / 2p r \/pw^ ' i^k dS = —ik / d t idS, 

J dK J dK 


where we have used Vpw^ = —ikd^pw^. The integral defining B(£,r) can be computed exactly, as shown in 
Remark n~75i below. 

Finally, G is the matrix of size (p,p) defined by 

G(£,m) = a K (pw m ,pwe), m,£ = l,...,p. 

The matrix G is Hermitian and invertible, provided that hxk is sufficiently small (see ©)• The matrix G can 
be expressed in terms of B and D as 

G = BD. 

Actually, G can be computed directly, without using quadrature formulae (see Remark 11.51 below). B and D 
have to be computed in any case. 

The matrix representation of the operator II is G~ l B , thus the matrix P of size (ukP^tikp) defined by 

P = DG~ 1 B 

is the matrix representation of the composition of the inclusion of V* K (K) in V Pk (K) after II. 
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1.5. The local PW-VEM bilinear form 

The local volume bilinear form a K (u, v ) cannot be computed exactly on general polygons K for u, v G V p (K). 
In a true PUM framework, once could consider as V(K) the space of generalized barycentric coordinates and 
use numerical quadrature. On the other hand, since the basis functions in V P (K) are highly oscillatory, the 
issue of the numerical integration is a non trivial one. 

Here we follow the VEM paradigm and we split the local bilinear form in a part that can be computed 
exactly (up to machine precision), and in a part that can be suitably approximated, provided that some stability 
properties are satisfied. Indeed, by using the projector n, the local bilinear form can be written as 

a K (u , v) = a K (Uu, Uv) + a K ((/ - U)u, (/ - U)v). 

The term a K (Uu, nu) can be computed in terms of the traces of the shape functions on dK only. Its matrix 
representation An is 

a u = b t g~ 1 b. 

On the contrary, the computation of the term a K ((/ — U)u, (/ — H)v) would require values of all the shape 
functions in the interior of K and it will be approximated (see Section [3] below). We denote its approximation 
by s K ((/ — n )u, (/ — n)u), and the local PW-VEM bilinear form can be written as 

a%{u, v) = a K (Uu, Uv) f s K ((/ - U)u, (/ - U)v). 

Since Hu* = u* for all u* G V*(K ), the following plane wave-consistency property holds true: 

= a K (u* p ,v) V«; G V;(K), V e V P (K). (17) 

1.6. The global PW-VEM formulation 

The global bilinear form defining the PW-VEM method is given by 


b h (u,v) = a h (u,v) 



uv d^S, 


(18) 


where 

a h (u,v)= a h(^,v)= [a K (Uu,Uv) + s K ((I-U)u,(I-n)v)], 

KeT h KeT h 

with 5 K (-, •) to be defined. Thus, the methods reads: find Uh p G V p (Th) such that 

bh(u hp ,v)=[ gvdS \/v G V p (Th)- (19) 

Jdn 

The boundary integral on the right-hand side of equation (H8|) is computed exactly (see Remark 11.51 below). 
thus the only integral which requires quadrature is the boundary integral on the right-had side of (fl9l) , containing 
the inhomogeneous boundary datum. In our theoretical analysis, nevertheless, in order to avoid complications, 
we assume that also this integral is computed exactly. 

Remark 1.4. We point out that the plane wave-consistency property m and the definition ofah( n -) imply 
that the Patch Test is satisfied, in the following sense. On any patch of elements, if the exact solution is a 
plane wave in one of the directions that define the local spaces V*(K) (or a linear combination of such plane 
waves), then the exact solution and the approximate solution coincide. 
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Remark 1.5. The computation of volume and edge integrals of products of plane waves, as well as that of edge 
integrals of products of plane waves by polynomial functions, can be done exactly (see J25 ]/ and f27\ pp. 20-21] 
In fact, if F is a mesh face (edge), denoting by a and b the coordinate vector of its endpoints, we have 

f e ifc(d m -dg)-x ^3 = | j p| e ifc(d m -d £ )-a f e ifc(d m -d*))-(b-a)t ^ _ 

J F JO 

= \F\e ik ^~ d ‘> a (ik( d m - de)) ■ (b - a)), 

where |F| denotes the length of F and, for z G C, 

_ , x [' zt , u {*—1. if z ^o 
®i(z)= e t dt = < 2 ? 

J ° y i if z = o. 

This formula also enters the computation of plane wave mass matrices, whose entries are integrals of the type 

/ e ifc(d m -d*)-x^y- 

Jk 

For m = I, this integral is simply \K\, the area of K. Whenever m ^ I, using Ae tkdx = —k 2 d ■de lkdyi , one 
has 


L 


p ik(dm-di) 


' dV = E 7 


FedK 


(d m — dj) • n F 
ik(d m - di) • (d m - djf) 


J e ik( 


where uf is the normal unit vector to F pointing outside K. 

For the matrix B, we have seen that B{I, r) = —ik f dK di • vk iprlFui dS. Thus, we need to compute integrals 
of the type 

j p J -(x)e ifc{d "*- d< >' x dS. 

If Vj is not an endpoint of F, then this integral is zero. Otherwise, denoting by a the coordinate vector of Vj 
and by b the coordinate vector of the other endpoint of F, we have 

f ^j(x)e iHdm - de) - x dS = \F\e ik(dm ~ de) ' a [ (1 - t)e ife(dm_d<)Hb_a)t di 

JF Jo 


= \F\e ik ^~ d ^ a $ 2 (ik(d m - d t )) • (b - a)), 


where 


z — 1 


1/2 


if z 7 ^ 0 
if z = 0 . 


$ 2 (z) = t (1 -t)e zt dt = 

Jo 

Similarly, for the integral on the right-hand side of CHI, we have to compute integrals of the type 

ipi(x)ipj ( X )e ife ( dm-d ^' x dS. 


L 


Whenever either Vi or Vj is not an endpoint of F, the integral is zero. Otherwise, we distinguish two cases. If 
i = j, denoting by a the coordinate vector of Vi = Vj and by b the coordinate vector of the other endpoint of F, 
we have 

[ (pi(-x)ipj(-x)e ik( ' dm ~ d ^' x dS = \F\e ik ( dm ~ de ' > ' a [ (1 - t ) 2 e *fe( d ™-d,))-(b-a) t dt 

Jf Jo 

= \F\e ik( ' dm ~ dt ' > ' a <S> 3 (ik(d m - d e )) • (b - a)), 
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where 

f i [ 2(e z - z-l)-zz , t n 

$ 3 (z)= / (l~t) 2 e zt dt= { zzz lfZ ^ 

J o [ 1/3 if z = 0. 

In the second case, when i ^ j, denoting by a and b the coordinate vector of Vj and Vi respectively, we have 


where 


f ipi(’x)ip j (x)e ik( ' dm ~ de ' > ' x dS = |i?|e ife(d -- de) ' a f (1 - t)te ik{ d -*- d *)Mb-a)t 
JF JO 

= \F\e ik ^ drn ~ dt J' a $ 4 (ik(d m - d,)) ■ (b - a)), 
$ 4 (s)= / (1 -t)te zi dt = 

Jo 


dt 


V(*-2) + s + 2 ifz ^ Q 


1/6 


if z = 0. 


2. Analysis 

The last step in defining our discretization scheme is the choice of the stabilization term s K ((I—U)u, (I—H)v). 
Before doing so, we first pose some abstract properties on the discrete problem that provide convergence results. 

2.1. Abstract result 

Let us introduce the /c-dependent norm for functions in 

ll< fcf n = l|V< n + fc 2 Ho, n , 

and the corresponding broken norm 

IMIl ,k,T h = IMIl ,k,K = (II^IU + IMIq.rM 

KeTh KeTh 


defined in the space iL 1 (7/ z ,) of broken iL 1 -functions. 

The continuous bilinear form b(-, •) satisfies the following continuity (see [571 Lemma 8.1.6]) and Garding 
inequality 


\b(u,v)\ < C CO nt |Mll,/c,0 Mil l,k,Q ’ 
Re[b(v, v)] + 2 k 2 ||-f ||q,o = \\ v \\i,k,n 


for all u, v G iL 1 (ll). 

Since for functions in the H-^ k Q -norm and the k ^-norm coincide, from here on, we will write 

(J* Hi k Th for both, whenever convenient. 

Theorem 2.1. Assume that the local stabilization forms s K ( •, •) are chosen in such a way that the following 
properties hold true: 

• continuity: there exists 7 > 0 such that, for all u,v G iL 1 (7/ l ); 

\a h {u,v)\ < 7l|w|li,fc,r h ll v lli,fe,r h : ( 20 ) 

• Garding inequality for the discrete operator; there exists a > 0 such that 

Re[b h (v, u)] + 2 k 2 ||w||o n > a ||v||i )fc)Th Vi; e V p (T h ). 


( 21 ) 
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Let u be the solution to problem m, and let Uh p be solution to the PW-VEM method (fT9|) in V p (Th)- Then, 
provided that h is small enough with respect to k (see ([32]) below), the following error estimate holds: 

1 + Q! + 7 / I, I, \ 

\\u-u h p\\ lkr <C - inf \\u-vi\\ lkT + inf L - , 

a \v I ev p (T h )" K p ev-(T h ) n np "i,k,T h I 


with C > 0 independent on h, k and p. Well-posedness of the PW-VEM method directly follows. 
Proof By the triangle inequality, we have 


||^ Uh p\\l,k,Th — W U Vl Wl,k,Th W Vl Uh pWl,k,Th 5 (^) 

for any vi G V p {Th )• We set Sh p = Uh p — vi, and proceed by estimating ||5^|| x k Th . 

The Garding inequality (|2T|) gives 

a W^ h p\\i,k,r h — R e [bh{Sh P ,Shpj} +2 k 2 ||^p||o ? Q =: I + IT (23) 

The term I can be treated as in |1J Th. 3.1]. Using the definition of 6/ l (-, •) and the discrete equation (fl9|) . 
for any v h P in v P 0h ), we get 


bh($h P ’> $h P ) — bh{ph p -, dh P ) O'hipi•> $h P ) ik I viSh P dS 

JdQ 

= / gSh P dS — aj£(vi,5h p ) — ik / vjSh p dS (insert v kp and use pw-consist. ([17]) ) 

J dn eTh, dn 

= f gShpdS- y] a%(vi - v* hp ,5 hp ) - a K (v* hp ,6 hp ) - ik f v{5 hp dS (use©) 
Jon KeTh KeT h 

= a{u, Shp) + ik / u Sh p d*S c^h{^i '^h P ’> ^h P ) ^ ^ & (Vh P i ^h P ) ik I viSh P d S 
J on k ^'T'h ” on 

= ^ a K (u-v* hp ,5 hp )+a h (v* hp -v I ,5 hp )+ik j (u-vVShpdS. 

KeT h ^ dn 


(24) 


In order to bound the last term, we observe that the trace inequality states that there exists a constant C > 0 
only depending on the shape of D such that, for all v G H 1 ( D), 


fcl/2 Nlo.an ^ Ckl/2 (diam(fi)- 1/2 |M| 0 n + |M|<^ ||Vn||J^) 
< C fc“ 1/2 diam(fi) -1 / 2 ||v|| ljfc>n + Q \\v\\ on + 
C (diam(0)“ 1/2 fc“ 1/2 |H| 1>fe>n + |M| 1)fcifl ) . 


< 


Due to the high-frequency assumption, it holds k 1 diam(D) 1 < 1/27T, and we conclude that there exists a 
constant C > 0 only depending on D such that, for all v G i7 1 (D), 


fcl/2 IMIo,an < C'NIi.jfe.n 
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Thus, since both {u — vi) and 5h p belong to we have 


ik 


[ (u-vi)S hp dS\ 
Jan 


<k\\u v i\\o^q ll^pllo,an 
— C \\u — vi 5 


(25) 


with a constant C > 0 only depending on 11. 

From ([24]) and ([25]) . using the continuity of the bilinear form and the continuity property ([20]) . we can 
conclude with the following estimate of the term /: 

I = R e[bh(8hp, $hp)] < Ci(l 7)(||r — v hp\\i^^ h + ||^ — v i\\i,k,T h ^ H^pIIi,/c,t^ • 

For the term //, we have 

11 = 2k \\S h p\\ 0 ^ < 2k ||£/ip|| 0) fi (11^ — u hp\\o,n + 11^ — v i llo,o) 

— ll^plli,fe,7^ (11^ — Uh p\\o,n + W u ~ llo,o) ’ 

due to k \\S hp \\ 0 n k Th . Inserting the bounds for / and II into ([23]) gives 

a IIMr,^ — Ci(l + 7)(|| w — v hp\\i^,T h ~ Vl Wi,k,Th) + 2/c(||n — Uhp\\ Q ,Q + 11 u ~ Fr|lo,o) 

< CII (1 'j) (\\u — v hp\\i^^ h + 11 ^ — v l\\l,k,Th) + \\u — UhpWo^Q 5 

where we have used again k ||u — u/|| 0 n < \\u — n/||i k j- . 

The term 11 u — Uh p || 0 n can be estimated by using a duality argument. Let us denote by ij; the solution of 
the dual problem 

b^Vy'i/j) = f v(u — Uh p )dV \/v G iL 1 (fl). (27) 

Jn 

Since 11 is assumed to be convex, ^ belongs to iL 2 (H) and satisfies 


IWIl,fc,Tfc - C W U Uh P Ho,0 ’ 

I1 2,0 — + k) || u — Uh p 11 0,0 ’ 


see m Prop. 8.1.4]. 

In correspondence to i/j G i7 2 (H), there exist G V*(Th) and ipi G V p (7h) and G V p (Th) such that 


V> “ ^hp\h,k,T h ~ C ( l + hk ^> h ll^lk k,Q > 
H-Mi, k ,r h <C(l + hk)h\m 2xn , 


with C > 0 independent of h, k and r ip. The first bound follows as in [2U Prop. 3.12 and 3.13]. The second one 
follows from combining [38l Th. 2.1] with the local approximation estimates in plane wave spaces. 
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= b(u — u hp , yj — ' 




Using (EH with v = u — Uh p , and inserting we have 

||u _ u hp\\l t n =H U - u hpA) = K u - u hp - V/) + - U hp ,ljj j) 

■^i) + / 

Jdn 

bh(u h p,^i) - / ^/dS' 

Jan 

= b(u - u hp , ^ - V/) - a(u hp , </>/) -ik u hp ^j dS 

Jan 

H ~ CLhiuhpi'ipi') ik / 'Uhp'ipi d^ 

Jan 

= 6(tx - u hp , ^ - Vh) + (ah(u hp , Vh) - a(u hp , Vh)) =: III + /V. 

The term ///is bounded using the continuity of the continuous form b(-, •), the second estimate in ([29]) and 
the regularity bounds ([28]) : 


III = b(u-u hp ^-^i) < C \\u — u hp \\ l k Th \\^-^i\\ 1AQ 
— Cm (1 + hk) h (1 + fc) ||r — 11^ — ^pllo,o * 

For the term /U, using the pw-consistency property (ED and the continuity of the discrete forms, we obtain 

/F = a h (u hp ,xlji) -‘a(u hp ,i/j I ) = Y (a%(u hp ,\l)i) ~ a K (u hp ,il)i)) 

KeT h 

= Y ( a h( u hp - v hp , V>l) - - <p, V’/)) (30) 

KeTh 

= Y ( a h( u hp - v* hp , Ipi - Vhp) - a K (u hp - v* hp , ipi - ip * hp )) 

KeT h 

< (l+i)\\u hp -v* hp \\^ Th ||^/H^p|| llfc ,r h - 

Next, using the bounds in ([29]) and ([28]) . we have 

< C (1 + hk) h (1 + fc) ||r - ^p||o,o J 

that we insert in ([30]) getting 

IV < Civ (1 + 7 ) (1 T hk) h (1 k) (\\u — u hp\\i^^ h + ||^ — v hp\\ 1 ^,%) ~ Uhp I1 0,0 * 

Therefore, we conclude with the following bound for \\u — Uh p \\ 0 

\\ u ~ u hp\\o^ C: (Cm + Civ) (1 + 7 ) (1 T- hk) h (1 V k) (\\u — u hp\\i^^ h + ||^ — Vh p\\i,k,Th^ (^) 

which, inserted into (126]) . gives 

a ll^p|l,fe,7fc ^Cii (1 + 7 )(\\u — v h P \\i^ : j- h + \\ u ~ v l\\l,k,Th) 

+ Cduai (1 + 7) k (1 + hk) h(l + k) (||r - Uh p \\ l k Th + \\u - v hp\\ l k Th ), 
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(Cduai = 2 (Cm + Civ)) and thus, owing to (|22j), 

<*\\u- Uhp\\ hktTh <Cn (l+a + 7)(||u-^ p || ljfejTfc + ||u~ 

+ ^duai (1 + 7) k (1 + hk) h(l + k) (||« - Uh p \\i t k,r h + lb _ ^plli.fc.rd 

Under the assumption 

C^duai (1 + 7) (1 + hk) hk (1 + k) < — , (32) 

we can take the \\u — Uh\h q- h term to the left-hand side and obtain 

W U - U hp\\ 1 , k ,T h — C “ {\\ U ~ V hp\\ ljkjTh + 11^ _ V l\\l,k,T h )i 

with C = 2C7/ + 1, which completes the proof. □ 

Remark 2.2. From ([32]) . it is clear that the threshold condition on h required in Theorem \2.1\ is that (1 + 
hk) hk (1 + k) be sufficiently small, which, in the relevant case of large k, is equivalent to requiring that hk 2 
be sufficiently small. This reflects the pollution effect of the h-version of the PW-VEM fffj. In fact, while a 
condition on hk is enough for the convergence of the best approximation (see Provo sition \2.S\ below), a stronger 
condition (namely, on hk 2 ) is required for the convergence of the method. 

The abstract convergence result of Theorem 12.11 combined with best approximation estimates of Helmholtz 
solutions within V p (Th) and V*(Th ), gives convergence rates. 

In order to state the following results, we define the weighted norm 

\\< k ,n = E k2(S ~ j) Kn, 

j=0 


where |-| . n denotes the standard seminorm in (Q). 

Proposition 2.3. Let u be a function in H i+1 (Ll), £ > 0, such that A u + k 2 u = 0 in Q. 
u hp ^ Vp(7h) and up G V p (Th), with p = 2m + 1, such that 

|| M - U hp\\l,k,T h - Cr l{hk) /l min{m ’° IMI min {m/}+l, fc.fi ’ 


ll u _ U /Ili,fc,r fc ^ Cr](hk)h n 


n{m,£} 


Ml 


min{ra,-£} + !,/c,Q ' 


Then there exist 


with C > 0 independent of h, k and u and 

r)(hk) = (1 + ( hk ) m+9 ) e(l~l p ) hk . 

Proof. The first bound follows from the local approximation estimates in plane wave spaces of [29] Th. 3.2.2] 
(see also j32J p. 831], [H]). The second bound can be obtained from the local approximation estimates in plane 
wave spaces by [38j Th. 2.1]. □ 

Theorem 12.11 and Proposition 12.31 immediately give the following convergence result. 

Corollary 2.4. Under the assumptions of Theorem \ 2. 11 if the solution u to problem (|2]) belongs to H^ +1 (Ul), 
£ > 1, then, provided that the threshold condition (EH) is satisfied, the following error estimate holds: 


\U - u hp \\ ltktTh < c V (hk) ||u|| minW}+1 ,*,n . 


with C > 0 independent on h, k and p = 2m + 1, and rj(hk) as in Proposition ^. 3\ 




16 


TITLE WILL BE SET BY THE PUBLISHER 


3. Stabilization term 


Let us give a sufficient condition on the stabilization term s K ((/ — II )u, (/ — II)u)) in order to guarantee the 
Garding inequality for the discrete operator. To this aim, we first state the following lemma. 

Lemma 3.1. For u G V P (K), we have 

k 2 [ (u-Uu)wdV = [ \7(u-Uu) - V^dV Vw e V*(K). 

Jk Jk 

Proof. The explicit form of the bilinear form that defines the projector II (see fl3J)) gives, for all w G V*(K), 


I 


VIRi • Vre dV 


— k 2 HuwdV = / Vu-'VwdV 
Jk Jk 


— k 2 uw dV, 
Jk 


□ 


that implies the assertion. 

Proposition 3.2. If the stabilization form satisfies the following condition: there exists a s > 0 such that, for 
all K G Th and v G V P (K), 

s K ((I - It)v, ( I - n)v) > a s ||V(J - n)u||o jjK ., (33) 

then the Garding inequality for the discrete operator holds true: 

Re[b h (v, u)] + 2 k 2 |Mlo n > min{o! s , 1} \\v\\\ XTh Vv G V p (T h ). 

Proof. We will make use of the trivial identity 

IK + V 2\\o,K = IKIIo,K + Iloilo, K + 2 ^ e [ v lV 2 dV . 

1 Jk j 


(34) 


Let us define a = min{o s , 1}. Due to a < 1, definition (fl8l) of b^(-, ■) and identity (f34|) with vy = Uv and 
V 2 = (/ — II)n, we have 

Re[Mu, v )\ + 2fc2 Iloilo,a = a h(v, v) + 2 k 2 |Mlo, a 

> Yi [aa K (Uv,Uv) + 2ak 2 \\nv\\ 2 0K 

Ken 1 


+ Y [s K (( I -R)v,(I-U)v)) + 2ak 2 \\(I-U)v\\l K 


Ker h 


V dak 2 Re [ [ Uv (I - U)v dR 
V Jk 


Ker h 


>aY [l|vnt;||^ + fc 2 ||n v ||2 J + Y Us\\V(I-H)v\\l Q + 2ak 2 \\(I-U)v\\l K 


KeTh 


Ker h 


V 2 aRe\ f \7Uv • V(/ - II)u dv] + V 2 ak 2 Re\ [ Uv(I-U) 
V JK J l JK 


Ker h 


Ker h 


)vdV 


where in the second inequality we have used ([33]) and Lemma [Til By using (f34|) with v\ — Vlln, v 2 = V(/—II)n, 
and then with v\ = Uv, v 2 = (I — II) u, we conclude 

Re[b h (v,v)] + 2k 2 \\v\\l n > min{o s ,l}(||Vn||o Q + k 2 \\v\\ 2 0 n ) = min{o s ,l} \\v\\l^ Th . 


□ 
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Under the assumption of Proposition 13.21 the Garding inequality of Theorem 12.11 holds true with a = 
min{<a s , 1 }. 

3.1. Choice ofs K (-,-) 

A first attempt in the choice of the stabilization term could be done by defining 

s K {{i -u)u,{i -u)v) = [ v(/-n)R- v(j-n)vdv, (35) 

JK 

so that (f33|) is satisfied with the equality and a s = 1. 

We notice that the continuity of the discrete bilinear form a/^.,.) follows from the continuity of the continuous 
bilinear forms and the continuity of the operator II (see Proposition HU). Therefore, flUD holds with 7 > 

1 + 2 (/ 3 mi 2 n + £min)> where Anin = VO.m. KeTh / 3 {h K k) = 

The matrix form of this stabilization term is 

(T^PfA(I - P), 

where where I is the identity matrix of appropriate size (number of directions p times number of mesh vertices 
%), P is the matrix representing the operator II defined in Section [L4l and A is the matrix of entries 

A(r,s)= [ VV’s-WrdV. 

JK 

However, when K is a generic polygon, we do not resort to explicit expressions of the basis functions ipj of 
V(K), and thus of ^ r , but we approximate (|35| ) . 

Adapting the rationale of VEM for elliptic problem to our case, we construct the (approximated) stabilization 
form as follows. Let ^ r (x) == pj (x) pwjg (x) and Vb( x ) = Pk ( x ) P w Km ( x ) be two basis functions in V P (K). Then 

VVv = (Vipj + ipj ik&e)pwji, 

= (V<p« + p K ik&m) pw^mi 


and thus 


* VVv = (Vp« • + ikp K d m • V(^j - ikifj di • - k 2 p K pj d m • d^)pw Km pw^. 


Taking into account the scaling of the terms in the brackets with respect to the elemental mesh size, we neglect 
the last three and replace the first by S K j/h 2 K . Therefore, we define s K ((/ — TP)u, (/ — n)p) in terms of its 
associated matrix Sk as 

S K =V Zr P) T M(I-P), (36) 

where M is the (scaled) plane wave mass matrix of size (ukP, n kP ) whose entries are M(r, s) = f K j^~ WJji dV . 
More precisely, if r = (j — l)p + £ and s = (k — l)p + m, 


M(r, s) 


xj) -ifed r (xx-x K ) 

n K 

/ pw m pw e dV 
Jk 

r 

° K j p ifcd m -(xjf-xj) -ifed r (xif-x ft ) 

/ i/c(d m -d £ )-(x- 

( 

( 

Jk 


( 37 ) 


We point out that the integral defining M(r, s) can be computed exactly; see Remark [131 
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On each element K , the elemental matrix of the complete PW-VEM is therefore 

B T G~ 1 B + JT^P) T M(I-P) + R (38) 

where D, B, G and P are defined in Section pi M is the scaled plane wave mass matrix defined in (l37|) , and 
R is the matrix associated with the bilinear form ik f dKnd q uv dS\ 

R(r,s)=ik / 'ip s 'ip r dS , r,s = 

JdKndn 

4. Numerical results 

Unless otherwise stated, in the following experiments we consider the problem © in the domain 11 = (0, l) 2 , 
and with boundary datum g such that the analytical solution is 

u(x) = H^\k |x - x 0 1), x 0 = (-0.25,0), 

where Hq 1 ^ is the zero-th order Hankel function of the first kind. In Figure [lj we report the real part of this 
solution u for the wave numbers k = 20, k = 40 and k = 60. All the errors reported in the following are relative 
errors in the L 2 -norm. 



Figure 1 . Real part of the function r(x) = H^^k |x — xo|), xo = (—0.25,0) in the domain 
H = (0, l) 2 , for k = 20 (left), k = 40 (center), and k = 60 (right). 


The function r(x) is a regular solution of the homogeneous Helmholtz problem, so that the h-convergence 
result of Corollary 12.41 provides order h m (with m — (p — l)/2) for the error in the H-l^ kl - h norm. An L 2 -error 
bound is then given by m- Because of the restriction (l32|) on the mesh size, we have hk 2 < C and the factor 
hk in ([3Tj) implies the gain of a factor h 1 / 2 for the L 2 -error, with respect to the ||.|| x k T ^-error. 

Due to the fact that high-order approximation properties of the plane wave spaces hold true only for solu¬ 
tions to the homogeneous Helmholtz equation, and our duality argument in the proof of Theorem 12.11 makes 
use of a dual problem with non zero right-hand side (see equation (l27]) h our analysis does not cover the asymp¬ 
totics in p. On the other hand, the p -version of the best approximation results (see [301 Th. 3.9, Rem. 3.13] 
and [HD Sect. 5.2] for the best approximation estimates in V*{Th ); for V p (Th ), use again [38l Th. 2.1]) give 
algebraic convergence (log (p)/p) £ , provided that the exact solution belongs to iL^ +1 (fl) and p is large enough, 
or exponential convergence, whenever the exact solution can be extended analytically outside D. Therefore, we 
also present numerical results for fixed h and increasing p. 
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First, we numerically evaluate in Section fl~Tl the loss of accuracy due to the approximation of the stabilization 
term. Then, in Section 14.21 we test the h- and p-versions of the PW-VEM on polygonal meshes, and compare 
the results obtained for different values of the wave number k. We conclude by testing the p-convergence for 
non analytic exact solutions (Section I4.2.3[) . 

4.1. Effects of the approximation of the stabilization term 

We recall that the choice s K ((/ — II)r, (/ — II)?;) = a K ((/ — II )u, (/ — II)?;) in the formulation (fl8l) coincides 
with the partition of unity (PUM) method, since the complete bilinear form is considered in this case. On 
triangular meshes, where the canonical VEM basis functions pj coincide with the classical Lagrangian Pi basis 
functions, we compare the error of the PUM, with that of the VEM with the stabilization (135|) (GRAD in the 
following) and that of the PW-VEM, whose elemental matrices are given by (l38|) . 

We have used structured triangular meshes, obtained from Cartesian subdivisions of U, dividing then each 
square into two triangles by one of the diagonals. We report in Table [1] the results obtained on four meshes 
containing, respectively, 8, 32, 128 and 512 triangles, using p = 13 plane waves per node, for the wave number 
k = 20. When the mesh size h is sufficiently small, GRAD is very close to PUM, showing that the sufficient 
requirement ([33]) on the stabilization term is indeed sharp. When the approximation of the stabilization term 
defined in (l36]) - (l37|) is used instead (PW-VEM), some accuracy is lost. However, the same order of convergence 
as for PUM is maintained. 


h 

PUM 

L 2 -error 

rate 

GRAI 

L 2 -error 

) 

rate 

PW-VE 

L 2 -error 

M 

rate 

7.0711e—01 

3.5355e—01 
1.7678e—01 
8.8388e—02 

1.9213e—02 

4.0683e—04 
3.4126e—06 
4.0164e—08 

5.5615 

6.8974 

6.4088 

7.1989e—01 
1.5517e—03 

3.3981e—06 
4.0148e—08 

8.8577 

8.8349 

6.4033 

4.1548e—01 

1.0990e—02 

1.2969e—04 
1.1089e—06 

5.2406 

6.4050 

6.8698 


Table 1. Relative error in the L 2 -norm for PUM, GRAD, PW-VEM for structured triangular 
meshes; k = 20, p = 13. 




Figure 2. Relative error in the L 2 -norm for PUM, GRAD, PW-VEM for the structured tri¬ 
angular mesh with 32 elements (h =3.5355e—01), for different values of p, for k = 20 (left) and 
km 40 (right). 

We also report in Figure [2] the errors with the three methods on the mesh with 32 elements, varying p, for 
k = 20 and k = 40, respectively. As in the h-error study, when the best approximation error is sufficiently 
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small (here, the number of plane waves per node is sufficiently large), GRAD coincides with PUM, and PW- 
VEM looses some accuracy, but still showing an exponential convergence behavior, with a comparable order of 
convergence. Notice that, when increasing the wave number fc, all these methods exhibit a larger preasymptotic 
region with slower convergence. Moreover, in both figures, and in many of the next plots, for large values of 
p, instability takes place due to the ill-conditioning of the plane wave basis, and the impact of roundoff error 
results in the increasing of the error; we refer to m for a similar behavior of the plane wave discontinuous 
Galerkin (PWDG) method. We remark that a similar phenomenon occurs for small h (see Section lT2l below i. 

We point out that the tests presented above do not provide a complete comparison between the PW-VEM 
and PUM (and/or GRAD), since we restricted to triangular meshes, where PUM can be exactly computed. 
A fair comparison should be carried out on more general meshes, where the use of PUM with generalized 
barycentric coordinates would need some quadrature formulas (we refer to m, where a PUM with generalized 
barycentric coordinates is compared to VEM with different stabilizations for the Poisson problem). On the 
other hand, these tests on triangular meshes suggest that there is a margin of improvement for the choice of 
the approximation of the stabilization term defined in (f36]) - ([37|) . 

4.2. Convergence 

We test convergence of the PW-VEM on (polygonal) Voronoi meshes made of 2 n elements, 3 < n < 9. We 
report in Figure 3 the meshes with 16 and 64 elements. 



Figure 3. Voronoi meshes with 16 (left) and 64 (right) elements. 



4.2.1. h-convergence 

We test first the h-convergence. We report in Figure (J4]) the relative errors in the L 2 -norm of the PW-VEM 
method on the Voronoi meshes, for pm 13 and k = 20, 40, 60. As pointed out at the beginning of Section 0J 
the expected convergence rate is m + 1/2, i.e., 6.5 for p = 13. This rate seems to be actually reached, after a 
pre-asymptotic region, which is wider for larger k. Also the results of Table [2] seem to confirm these rates (6.5 
for p — 13 and 7.5 for p = 15). However, since the sequence of Voronoi meshes considered here is not made by 
nested meshes, the maximum of the diameters of the elements K G Th (definition of h for Table [2]) does not 
vary by a factor from one mesh to the other, so that a precise assessment of the convergence rate is difficult. 
Nevertheless, we can conclude that high order rate is provided also on polygonal meshes. We also report in 
Figure [5] the error plots for k = 20 and for different values of p. For fine meshes and high p, instability takes 
place and the error increases, as already observed. 

In order to test the pollution effect, we have run the PW-VEM, with p = 9, on the Voronoi meshes made of 2 n 
elements, 3 < n < 9, and k chosen such that the product hk is the constant 3 (i.e., k = 6.26, 8.70, 13.63, 19.51, 
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Figure 4. Relative error in the L 2 -norm of the PW-VEM for p = 13 and k = 20, 40, 60 on 
the Voronoi meshes. 


h 

p = i; 
L 2 -error 

3 

rate 

p = 11 

L 2 -error 

rate 

3.4487e—01 
2.2017e—01 
1.5376e—01 
1.1806e—01 

8.4571e—02 

2.7882e—02 
4.6014e—03 
4.0962e—04 
4.1264e—05 
4.1597e—06 

4.0147 

6.7374 

8.6874 

6.8785 

1.1374e—02 

1.5253e—03 

6.6821e—05 

5.3076e-06 
9.3361e—07 

4.4772 

8.7123 

9.5868 

5.2096 


Table 2. Relative error in the L 2 -norm for PW-VEM for k = 20, with p — 13 and p = 15, on 
Voronoi meshes. 


25.41, 35.47, 51.59, respectively). The results reported in Figure [6] confirm that the h -version of the PW-VEM 
is affected by pollution, i.e., the boundedness of the product hk is not sufficient to guarantee convergence of the 
discretization error. 

4.2.2. p-convergence 

In Figure [3 we report the relative errors in the L 2 -norm for different values of p in the case k = 20, on the 
Voronoi meshes with 16 and 64 elements. We see exponential convergence in p, before instability for high p 
takes place, as in Figure [2j left, for the triangular meshes. In the case of larger elements (Figure [3 left), the 
pre-asymptotic region is wider, while for smaller elements (Figure [3 right), the instability starts to take place 
for smaller values of p. Notice that the considered meshes also contain edges which are small, as compared to 
the mesh size. These results, together with those in Section 14.2.11 thus suggest that the PW-VEM is robust in 
case of degenerating sides. 

In Figure[8j a polygonal mesh with non convex elements (100 elements) is represented, and the corresponding 
L 2 -error for k = 20 and different values of p is plotted. The results are comparable to those obtained with meshes 
with convex elements, confirming robustness of the method also in the presence of non convex elements. Due 
to the presence of small elements, instability starts to take place for smaller values of p, as compared to the 
previous, more regular, meshes. 
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Figure 5. Relative error in the L 2 -norm of the PW-VEM for k = 20 and different values of p 
on the Voronoi meshes. 



Figure 6 . Relative error in the L 2 -norm of the PW-VEM on the Voronoi meshes for hk = 3 
and for p = 9. 


We compare now the convergence behaviour of the PW-VEM for different values of the wavenumber k. We 
plot in Figure [9] the L 2 -error of the PW-VEM on the Voronoi mesh with 64 element, for different values of p 
and for k = 20, 40, 60. For larger values of fc, the pre-asymptotic region is wider, while instability starts to take 
place for larger values of p. 

4.2.3. p-convergence for non smooth solutions 

Finally, we test the p-convergence in the case of solutions with low regularity. To this aim, as in [30, Sect 4], 
we consider now the problem 0 again in the domain Q = (0, l) 2 , and with boundary datum g such that the 
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Figure 7. Relative error in the L 2 -norm for PW-VEM for k = 20 on Voronoi meshes with 16 
(left) and 64 (right) elements, for different values of p. 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 8. Polygonal mesh with 100 elements containing non convex elements (left), and rel¬ 
ative error in the L 2 -norm for PW-VEM for k = 20 and different values of p (right). 


analytical solution is given, in polar coordinates x = (r cos'd, r sin$), by 

r(x — xo) = J ^( kr ) cos(£$), £ > 0, xo = (0,0.5), (39) 

where J £ is the Bessel function of the first kind and order £. We choose k = 10. For integer £, u can be extended 
analytically outside fi, otherwise its derivatives have a singularity at xo- We report in Figure [TO] the real part 
of u for £ = 1 (regular case), £ = 3/2 (u E H 2 (Q)) and £ = 2/3 (u 0 H 2 (Q)). 

We report in Figure [IT] the relative errors in the L 2 -norm for the three cases, obtained on the structured 
triangular mesh with 32 elements, for different values of p. While for the smooth solution (£ = 1) we have 
exponential convergence, before instability takes place, as in the previous experiments, in the other two cases 
the convergence is algebraic, even if the rates are unclear, the results for £ = 3/2 being better than those for 
£ = 2/3. Similar results were obtained in [30] for the PWDG method. We also report the results obtained with 
the PUM, for the sake of comparison (this is the reason why we have used triangular meshes). 
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Figure 9. Relative error in the L 2 -norm of the PW-VEM on the Voronoi mesh with 64 ele¬ 
ments, for different values of p and for k = 20, 40, 60. 



0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 


Figure 10. Real part of the function r(x — xo) = J^(kr) cost'd), xo = (0,0.5) in the domain 
Q = (0, l) 2 , for k = 10 and £ = 1 (left), £ = 3/2 (center), and £ = 2/3 (right). 

5. Conclusions 

We have presented a first PW-VEM scheme for the discretization of the Helmholtz equation. H 1 -conformity 
is guaranteed by the VEM framework, while high order convergence for the homogeneous problem, for smooth 
analytical solutions, is achieved through a plane wave enrichment of the approximating spaces. An h- version 
error analysis of the PW-VEM is derived, while numerical results also show exponential convergence of the 
p- version for smooth analytical solutions. These preliminary results highlight that a suitable interplay of h, p 
and k is important in order to obtain quasi-optimality, and suggest that similar results as those of the complete 
hp- analysis of SOI could hold. 

Several restrictions on the model problem we have considered could be easily removed, and the method could 
be extended from 2D to 3D, and to acoustic scattering problems. The extension of the method and its analysis 
to problems with non constant coefficients or non zero source terms is also an interesting issue. Exploring 
alternative choices for the stabilization term, the projection operator n and/or of the space onto which to 
project, as well as extensions to non uniform p and higher order VEM will be subject of future research. 
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p/ log (p) 

Figure 11. Relative error in the L 2 -norm of the PW-VEM and PUM on the structured tri¬ 
angular mesh with 32 elements, for the model problem with exact solution ([39]) for k = 10 and 
£ = 1 (top, left), £ = 3/2 (top, right), and £ = 2/3 (bottom), for different values of p. 
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